Legendre et Gallagher 2001

Create data of Fig. 3a

fig3a_data <- data.frame(
  site = seq(1:19),
  sp1 = c(7,4,2,1, rep(0,15)),
  sp2 = c(1,2,4,7,8,7,4,2,1, rep(0,10)),
  sp3 = c(rep(0,5),1,2,4,7,8,7,4,2,1,rep(0,5)),
  sp4 = c(rep(0,10),1,2,4,7,8,7,4,2,1),
  sp5 = c(rep(0,15), 1,2,4,7),
  sp6 = c(0,1,1,rep(0,16)),
  sp7 = c(rep(0,6),2,1,rep(0,11)),
  sp8 = c(rep(0,11),3,1,rep(0,6)),
  sp9 = c(rep(0,16), 4,1,0)
)

pos <- position_dodge(0.3)

fig3a_data %>%
  gather(species, abundance, -site) %>%
  ggplot(data= ., aes(x = site, y = abundance, fill = species, shape = species, group = species)) +
  geom_line(aes(col = species), position = pos, alpha = 0.6) +
  geom_point(data = . %>% filter(abundance != 0), alpha=0.6, size = 2.5, position = pos) +
  scale_shape_manual(values=rep(c(21:25),4)) +
  scale_x_continuous(breaks = seq(1:19)) +
  theme_linedraw()

Figure 1:

fig3a_data <- fig3a_data %>%
  column_to_rownames("site")

Reproduce Fig. 4

Using decostand :

• profile transformation: Y.tr = decostand(Y, “total”)
• chord transformation: Y.tr = decostand(Y, “norm”)
• log-chord transformation: Y.tr = decostand(log1p(Y), “norm”)
• Hellinger transformation: Y.tr = decostand(Y, “hellinger”)
• chi-square transformation: Y.tr = decostand(Y, “chi.sq”)

dist.ldc can also be used

comm <- fig3a_data

# Figure 4a : PCA on raw data
#-----------------------------

ord <- rda(comm)
fig_4a <- plot_ord(ord, comm) +
  ggtitle("PCA on raw data")

# Figure 4b : CA on raw data
#---------------------------

ord <- cca(comm)
fig_4b <- plot_ord(ord, comm) +
  ggtitle("CA on raw data")

# Figure 4c : PCA on chi2-transformed data
#-----------------------------------------

comm_tr <- decostand(comm,"chi.sq")
ord <- rda(comm_tr)
fig_4c <- plot_ord(ord, comm) +
  ggtitle("PCA on chi2-transformed data")

# Figure 4d : PCA on chord-transformed data
#------------------------------------------

comm_tr <- decostand(comm,"norm")
ord <- rda(comm_tr)
fig_4d <- plot_ord(ord, comm) +
  ggtitle("PCA on chord-transformed data")

# Figure 4e : PCA on species-profiles transformed data
#-----------------------------------------------------

comm_tr <- decostand(comm,"total")
ord <- rda(comm_tr)
fig_4e <- plot_ord(ord, comm) +
  ggtitle("PCA on species-profiles transformed data")

# Figure 4f : PCA on Hellinger transformed data
#---------------------------------------------

comm_tr <- decostand(comm,"hellinger")
ord <- rda(comm_tr)
fig_4f <- plot_ord(ord, comm) +
  ggtitle("PCA on Hellinger-transformed data")

# Figure 4g : PCoA with bray-curtis distance
#---------------------------------------------

comm_d <- vegdist(comm, "bray")
ord <- cmdscale(comm_d)
fig_4g <- plot_ord(ord, comm, reverse_axis2 = TRUE) +
  ggtitle("PCoA with Bray-curtis")

(fig_4a + fig_4b + fig_4c) / 
  (fig_4d + fig_4e + fig_4f) / 
  (fig_4g + plot_spacer() + plot_spacer()) +
  plot_annotation(tag_levels = 'a', tag_prefix="(", tag_suffix = ")")

Figure 2:

• Add variance explained to axes labels
• Do log-chord

Comparison of Fig. 4 with JSDM

Comparison with boral

Normal

# Normal - without row effects
#----------------------------

boral_normal <- boral(y = comm, family = "normal", lv.control = list(num.lv = 2), row.eff = "none")

boral_p1 <- plot_boral(boral_normal)+
  scale_y_reverse() +
  scale_x_reverse() +
  ggtitle("Normal - without row effects")

# Normal - fixed row effects
#----------------------------

boral_normal_fixed <- boral(y = comm, family = "normal", lv.control = list(num.lv = 2), row.eff = "fixed")

boral_p2 <- plot_boral(boral_normal_fixed)+
  scale_y_reverse() +
  ggtitle("Normal - fixed row effects")

# Normal - random row effects
#----------------------------

boral_normal_rand <- boral(y = comm, family = "normal", lv.control = list(num.lv = 2), row.eff = "random")

boral_p3 <- plot_boral(boral_normal_rand)+
  scale_y_reverse() +
  scale_x_reverse() +
  ggtitle("Normal - random row effects")

# Hellinger - Normal - without row effects
#----------------------------
comm_tr <- decostand(comm,"hellinger")

boral_normal_hell <- boral(y = comm_tr, family = "normal", lv.control = list(num.lv = 2), row.eff = "none")

boral_p4 <- plot_boral(boral_normal_hell)+
  scale_y_reverse() +
  scale_x_reverse() +
  ggtitle("Hellinger-transformed - Normal - without row effects")

# Normal - fixed row effects
#----------------------------

boral_normal_hell_fixed <- boral(y = comm_tr, family = "normal", lv.control = list(num.lv = 2), row.eff = "fixed")

boral_p5 <- plot_boral(boral_normal_hell_fixed)+
  scale_y_reverse() +
  ggtitle("Hellinger-transformed - Normal - fixed row effects")

# normal - random row effects
#----------------------------

boral_normal_hell_rand <- boral(y = comm_tr, family = "normal", lv.control = list(num.lv = 2), row.eff = "random")

boral_p6 <- plot_boral(boral_normal_hell_rand)+
  scale_y_reverse() +
  scale_x_reverse() +
  ggtitle("Hellinger-transformed - Normal - random row effects")

(boral_p1 + boral_p4) / (boral_p2 + boral_p5) / (boral_p3 + boral_p6) +
  plot_annotation(tag_levels = 'a', tag_prefix="(", tag_suffix = ")")

Figure 3:

(boral_p1 + boral_p2 + boral_p3)+
  plot_annotation(tag_levels = 'a', tag_prefix="(", tag_suffix = ")")

Figure 4:

Poisson & Negative binomial

# Poisson - without row effects
#----------------------------

boral_poisson <- boral(y = comm, family = "poisson", lv.control = list(num.lv = 2), row.eff = "none")

boral_p1 <- plot_boral(boral_poisson)+
  scale_y_reverse() +
  scale_x_reverse() +
  ggtitle("Poisson - without row effects")

# Poisson - fixed row effects
#----------------------------

boral_poisson_fixed <- boral(y = comm, family = "poisson", lv.control = list(num.lv = 2), row.eff = "fixed")

boral_p2 <- plot_boral(boral_poisson_fixed)+
  scale_y_reverse() +
  ggtitle("Poisson - fixed row effects")

# Poisson - random row effects
#----------------------------

boral_poisson_rand <- boral(y = comm, family = "poisson", lv.control = list(num.lv = 2), row.eff = "random")

boral_p3 <- plot_boral(boral_poisson_rand)+
  scale_y_reverse() +
  scale_x_reverse() +
  ggtitle("Poisson - random row effects")

# Negative binomial - without row effects
#--------------------------------------

boral_negbin <- boral(y = comm, family = "negative.binomial", lv.control = list(num.lv = 2), row.eff = "none")

boral_p4 <- plot_boral(boral_negbin) +
  scale_y_reverse() +
  ggtitle("Negative binomial - without row effects")

# Negative binomial - fixed row effects
#--------------------------------------

boral_negbin_fixed <- boral(y = comm, family = "negative.binomial", lv.control = list(num.lv = 2), row.eff = "fixed")

boral_p5 <- plot_boral(boral_negbin_fixed) +
  scale_y_reverse() +
  ggtitle("Negative binomial - fixed row effects")

# Negative binomial - random row effects
#--------------------------------------

boral_negbin_rand<- boral(y = comm, family = "negative.binomial", lv.control = list(num.lv = 2), row.eff = "random")

boral_p6 <- plot_boral(boral_negbin_rand) +
  scale_y_reverse() +
  ggtitle("Negative binomial - random row effects")

(boral_p1 + boral_p4) / (boral_p2 + boral_p5) / (boral_p3 + boral_p6) +
  plot_annotation(tag_levels = 'a', tag_prefix="(", tag_suffix = ")")

Figure 5:

Residuals of the models

• Normal - without row effects

plot(boral_normal, ask = F, mfrow = c(2,2))

## NULL

Figure 6:

Figure 7:

Figure 8:

Figure 9:

• Poisson - without row effects

plot(boral_poisson, ask = F, mfrow = c(2,2))

## NULL

Figure 10:

Figure 11:

Figure 12:

Figure 13:

• Negative binomial - without row effects

plot(boral_negbin, ask = F, mfrow = c(2,2))

## NULL

Figure 14:

Figure 15:

Figure 16:

Figure 17:

Comparison with mvabund

meanvar.plot(comm)

Figure 18:

Session info

sessionInfo()

## R version 3.5.0 (2018-04-23)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS  10.15.5
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] boral_1.8       coda_0.19-3     mvabund_4.0.1   vegan_2.5-6    
##  [5] lattice_0.20-35 permute_0.9-5   ggboral_0.1.7   patchwork_0.0.1
##  [9] ggrepel_0.8.2   forcats_0.3.0   stringr_1.4.0   dplyr_1.0.2    
## [13] purrr_0.3.2     readr_1.1.1     tidyr_1.0.0     tibble_2.1.3   
## [17] ggplot2_3.3.0   tidyverse_1.2.1
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-137     lubridate_1.7.4  httr_1.4.1       tools_3.5.0     
##  [5] backports_1.1.4  R6_2.4.1         R2WinBUGS_2.1-21 mgcv_1.8-23     
##  [9] questionr_0.7.0  colorspace_1.4-1 withr_2.1.2      tidyselect_1.1.0
## [13] compiler_3.5.0   cli_2.0.2        rvest_0.3.4      xml2_1.2.2      
## [17] labeling_0.3     bookdown_0.13    scales_1.0.0     mvtnorm_1.1-0   
## [21] digest_0.6.25    rmarkdown_1.15   pkgconfig_2.0.3  htmltools_0.3.6 
## [25] highr_0.7        rlang_0.4.8      readxl_1.1.0     rstudioapi_0.9.0
## [29] shiny_1.3.2      generics_0.0.2   jsonlite_1.6     magrittr_1.5    
## [33] Matrix_1.2-14    Rcpp_1.0.5       munsell_0.5.0    fansi_0.4.1     
## [37] abind_1.4-5      lifecycle_0.2.0  stringi_1.4.6    yaml_2.2.0      
## [41] MASS_7.3-51.4    plyr_1.8.4       fishMod_0.29     grid_3.5.0      
## [45] parallel_3.5.0   promises_1.0.1   crayon_1.3.4     miniUI_0.1.1.1  
## [49] haven_2.1.1      hms_0.4.2        knitr_1.24       pillar_1.4.2    
## [53] boot_1.3-20      corpcor_1.6.9    reshape2_1.4.3   codetools_0.2-15
## [57] R2jags_0.5-7     glue_1.4.2       evaluate_0.14    modelr_0.1.5    
## [61] vctrs_0.3.5      rmdformats_0.3.5 httpuv_1.5.2     cellranger_1.1.0
## [65] gtable_0.3.0     assertthat_0.2.1 xfun_0.9         mime_0.7        
## [69] tweedie_2.3.2    xtable_1.8-4     broom_0.5.2      later_0.8.0     
## [73] rjags_4-10       cluster_2.0.7-1  statmod_1.4.34   ellipsis_0.3.0